Setup

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(clusterProfiler)
## 
## clusterProfiler v4.10.1  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
## 
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
## 
## Attaching package: 'clusterProfiler'
## 
## The following object is masked from 'package:purrr':
## 
##     simplify
## 
## The following object is masked from 'package:stats':
## 
##     filter
library(org.Hs.eg.db)
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min
## 
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Loading required package: IRanges
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     rename
## 
## The following objects are masked from 'package:lubridate':
## 
##     second, second<-
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following object is masked from 'package:utils':
## 
##     findMatches
## 
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     slice
## 
## The following object is masked from 'package:lubridate':
## 
##     %within%
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## The following object is masked from 'package:grDevices':
## 
##     windows
## 
## 
## Attaching package: 'AnnotationDbi'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     select
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(enrichplot)
library(pathview)
## ##############################################################################
## Pathview is an open source software package distributed under GNU General
## Public License version 3 (GPLv3). Details of GPLv3 is available at
## http://www.gnu.org/licenses/gpl-3.0.html. Particullary, users are required to
## formally cite the original Pathview paper (not just mention it) in publications
## or products. For details, do citation("pathview") within R.
## 
## The pathview downloads and uses KEGG data. Non-academic uses may require a KEGG
## license agreement (details at http://www.kegg.jp/kegg/legal.html).
## ##############################################################################
library(DOSE)
## DOSE v3.28.2  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
## 
## If you use DOSE in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609
knitr::opts_chunk$set(fig.width = 10, dpi = 300)

dir.create(path = "./results/FA", showWarnings = FALSE)

Results

TLR10 light vs. dark 3h

Gene ontology

out.dir <- "./results/FA/TLR10_light_3h_vs_TLR10_dark_3h/"

dir.create(out.dir)

res.TLR3h <- read.csv("results/DESeq/TLR10_light_3h_vs_TLR10_dark_3h_unfiltered.csv") %>%
  as_tibble() %>%
  arrange(desc(log2FoldChange))

gene_list.TLR3h <- res.TLR3h$log2FoldChange
names(gene_list.TLR3h) <- res.TLR3h$ENSG
  
set.seed(5892)
gseGO.TLR3h <- gseGO(gene_list.TLR3h,
                     OrgDb = org.Hs.eg.db,
                     keyType = "ENSEMBL",
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = TRUE)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseGO.TLR3h <- setReadable(gseGO.TLR3h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLR3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      BP 
## #...@keytype      ENSEMBL 
## #...@geneList     Named num [1:32287] 1.61 1.54 1.4 1.37 1.31 ...
##  - attr(*, "names")= chr [1:32287] "ENSG00000204388" "ENSG00000204389" "ENSG00000051108" "ENSG00000175197" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...49 enriched terms found
## 'data.frame':    49 obs. of  11 variables:
##  $ ID             : chr  "GO:0006457" "GO:0051604" "GO:0002181" "GO:0061077" ...
##  $ Description    : chr  "protein folding" "protein maturation" "cytoplasmic translation" "chaperone-mediated protein folding" ...
##  $ setSize        : int  200 495 154 71 134 167 156 41 121 98 ...
##  $ enrichmentScore: num  0.676 0.525 0.657 0.769 0.671 ...
##  $ NES            : num  2.01 1.66 1.9 2.02 1.92 ...
##  $ pvalue         : num  2.00e-09 8.10e-08 5.01e-07 1.27e-06 2.65e-06 ...
##  $ p.adjust       : num  9.65e-06 1.95e-04 8.05e-04 1.52e-03 2.55e-03 ...
##  $ qvalue         : num  8.33e-06 1.68e-04 6.95e-04 1.32e-03 2.20e-03 ...
##  $ rank           : num  3165 4282 3050 2216 3646 ...
##  $ leading_edge   : chr  "tags=44%, list=10%, signal=40%" "tags=35%, list=13%, signal=31%" "tags=60%, list=9%, signal=55%" "tags=44%, list=7%, signal=41%" ...
##  $ core_enrichment: chr  "HSPA1B/HSPA1A/HSPA1L/DNAJB1/HSPB1/HSPH1/PPID/DNAJA1/GRPEL2/CHORDC1/ENGASE/DNAJB2/HSPA5/ARL2/AIP/LMAN2L/ST13/CDC"| __truncated__ "HSPA1B/HSPA1A/HSPA1L/DNAJB1/CHAC1/HSPB1/HSPH1/SIRT4/PPID/DNAJA1/GRPEL2/CHORDC1/ENGASE/DNAJB2/HSPA5/ARL2/AIP/NOL"| __truncated__ "EIF4A2/RPL28/RPS19/AARS1/UBA52/EIF3H/RPL30/RPL13A/RPS12/RPS4X/RPL31/RPLP1/RPS28/RPS15A/RPL37A/RPS9/RPL14/RPL13/"| __truncated__ "HSPA1B/HSPA1A/HSPA1L/DNAJB1/HSPB1/HSPH1/PPID/CHORDC1/DNAJB2/HSPA5/ST13/SDF2L1/SDF2/FKBP2/PFDN5/PPIB/FKBP4/PFDN4"| __truncated__ ...
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseGO.TLR3h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)

TLR3h.go.plot <- nrow(gseGO.TLR3h) > 0
dotplot(gseGO.TLR3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseGO.TLR3h, showCategory = 5, color.params = list(foldChange = gene_list.TLR3h),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

gseGO.TLR3h.pwts <- pairwise_termsim(gseGO.TLR3h)

emapplot(gseGO.TLR3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

KEGG (Kyoto Encyclopedia of Genes and Genomes)

ids <- bitr(names(gene_list.TLR3h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
## 'select()' returned 1:many mapping between keys and columns
ids <- ids[!duplicated(ids$ENSEMBL),]

gene_list.TLR3h.entrez <- gene_list.TLR3h[names(gene_list.TLR3h) %in% ids$ENSEMBL]
names(gene_list.TLR3h.entrez) <- ids$ENTREZID
gene_list.TLR3h.entrez <- gene_list.TLR3h.entrez[!duplicated(names(gene_list.TLR3h.entrez))]
set.seed(5923)
gseKEGG.TLR3h <- gseKEGG(gene_list.TLR3h.entrez,
                         organism = "hsa",
                         keyType = "ncbi-geneid",
                         eps = 0,
                         maxGSSize = 500,
                         minGSSize = 15,
                         seed = T)
## Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"...
## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/hsa"...
## Reading KEGG annotation online: "https://rest.kegg.jp/conv/ncbi-geneid/hsa"...
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseKEGG.TLR3h <- setReadable(gseKEGG.TLR3h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLR3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     hsa 
## #...@setType      KEGG 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 1.61 1.54 1.4 1.37 1.31 ...
##  - attr(*, "names")= chr [1:22414] "3304" "3303" "9709" "1649" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...73 enriched terms found
## 'data.frame':    73 obs. of  11 variables:
##  $ ID             : chr  "hsa03010" "hsa05171" "hsa04141" "hsa05020" ...
##  $ Description    : chr  "Ribosome" "Coronavirus disease - COVID-19" "Protein processing in endoplasmic reticulum" "Prion disease" ...
##  $ setSize        : int  129 203 164 255 57 131 255 149 54 292 ...
##  $ enrichmentScore: num  0.786 0.597 0.63 0.567 0.791 ...
##  $ NES            : num  2.49 1.98 2.04 1.93 2.22 ...
##  $ pvalue         : num  8.61e-17 3.21e-09 6.93e-09 6.23e-09 5.65e-08 ...
##  $ p.adjust       : num  2.87e-14 5.37e-07 5.79e-07 5.79e-07 3.77e-06 ...
##  $ qvalue         : num  1.87e-14 3.48e-07 3.76e-07 3.76e-07 2.45e-06 ...
##  $ rank           : num  1880 2592 2319 2411 2044 ...
##  $ leading_edge   : chr  "tags=62%, list=8%, signal=57%" "tags=41%, list=12%, signal=37%" "tags=36%, list=10%, signal=32%" "tags=34%, list=11%, signal=30%" ...
##  $ core_enrichment: chr  "RPL28/RPS19/UBA52/RPL30/MRPL34/RPL13A/RPS12/RPS4X/RPL31/RPLP1/RPS28/RPS15A/RPL37A/RPS9/RPL14/RPL13/RPL37/RPS11/"| __truncated__ "RPL28/RPS19/UBA52/RPL30/RPL13A/RPS12/RPS4X/RPL31/RPLP1/RPS28/RPS15A/RPL37A/RPS9/RPL14/RPL13/RPL37/RPS11/RSL24D1"| __truncated__ "HSPA1B/HSPA1A/HERPUD1/DDIT3/HSPA1L/DNAJB1/ATF4/HSPH1/XBP1/NFE2L2/DNAJA1/DNAJB2/HSPA5/TMEM258/SSR4/UBE2D1/LMAN2/"| __truncated__ "HSPA1B/HSPA1A/DDIT3/HSPA1L/ATF4/HSPA5/NDUFA3/COX5B/COX4I1/COX7C/UQCR11/UQCRB/NDUFS8/ATF2/ATP5F1E/NDUFA2/UQCRQ/C"| __truncated__ ...
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseKEGG.TLR3h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)

TLR3h.kegg.plot <- nrow(gseKEGG.TLR3h) > 0
dotplot(gseKEGG.TLR3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseKEGG.TLR3h, showCategory = 5, color.params = list(foldChange = gene_list.TLR3h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseKEGG.TLR3h.pwts <- pairwise_termsim(gseKEGG.TLR3h)

emapplot(gseKEGG.TLR3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

pathview(gene.data = gene_list.TLR3h.entrez, pathway.id = gseKEGG.TLR3h@result$ID[1], species = "hsa")
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq
## Info: Writing image file hsa03010.pathview.png
filename <- paste0(gseKEGG.TLR3h@result$ID[1], ".pathview.png")

file.rename(filename, paste0(out.dir, filename))
## [1] TRUE
knitr::include_graphics(paste0(out.dir, filename))

Disease ontology

set.seed(5534)
gseDO.TLR3h <- gseDO(gene_list.TLR3h.entrez,
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseDO.TLR3h <- setReadable(gseDO.TLR3h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLR3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      DO 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 1.61 1.54 1.4 1.37 1.31 ...
##  - attr(*, "names")= chr [1:22414] "3304" "3303" "9709" "1649" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...21 enriched terms found
## 'data.frame':    21 obs. of  11 variables:
##  $ ID             : chr  "DOID:1470" "DOID:2841" "DOID:12306" "DOID:1176" ...
##  $ Description    : chr  "major depressive disorder" "asthma" "vitiligo" "bronchial disease" ...
##  $ setSize        : int  37 255 70 264 19 17 111 24 158 300 ...
##  $ enrichmentScore: num  0.813 0.512 0.693 0.499 0.902 ...
##  $ NES            : num  2.11 1.73 2.01 1.7 2.08 ...
##  $ pvalue         : num  3.82e-06 4.72e-06 1.01e-05 2.23e-05 4.03e-05 ...
##  $ p.adjust       : num  0.00226 0.00226 0.00323 0.00534 0.00594 ...
##  $ qvalue         : num  0.00205 0.00205 0.00293 0.00485 0.00538 ...
##  $ rank           : num  1299 1977 1359 1977 583 ...
##  $ leading_edge   : chr  "tags=14%, list=6%, signal=13%" "tags=13%, list=9%, signal=12%" "tags=23%, list=6%, signal=22%" "tags=12%, list=9%, signal=12%" ...
##  $ core_enrichment: chr  "HSPA1A/HSPA1L/COMT/PTEN/EIF4B" "HSPA1B/HSPA1A/HSPB1/NFE2L2/LTA4H/UCN/PRMT2/CAT/GSTP1/SOD1/GSTM3/SELPLG/IFNGR1/MIF/ORMDL1/HLA-C/ARG2/NQO1/ANG/LG"| __truncated__ "HSPA1A/XBP1/NFE2L2/CAT/COMT/GSTP1/TSLP/HLA-A/PSMB8/GPX1/HLA-C/PSMB9/MSRB2/HLA-B/FAS/PTGS2" "HSPA1B/HSPA1A/HSPB1/NFE2L2/LTA4H/UCN/PRMT2/CAT/GSTP1/SOD1/GSTM3/SELPLG/IFNGR1/MIF/ORMDL1/HLA-C/ARG2/NQO1/ANG/LG"| __truncated__ ...
## #...Citation
##   Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
##   R/Bioconductor package for Disease Ontology Semantic and Enrichment
##   analysis. Bioinformatics 2015, 31(4):608-609
write.csv(data.frame(gseDO.TLR3h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)

TLR3h.DO.plot <- nrow(gseDO.TLR3h) > 0
dotplot(gseDO.TLR3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseDO.TLR3h, showCategory = 5, layout = "dh", color.params = list(foldChange = gene_list.TLR3h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseDO.TLR3h.pwts <- pairwise_termsim(gseDO.TLR3h)

emapplot(gseDO.TLR3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

TLR10 light vs. dark 16h

Gene ontology

out.dir <- "./results/FA/TLR10_light_16h_vs_TLR10_dark_16h/"

dir.create(out.dir)

res.TLR16h <- read.csv("results/DESeq/TLR10_light_16h_vs_TLR10_dark_16h_unfiltered.csv") %>%
  as_tibble() %>%
  arrange(desc(log2FoldChange))

gene_list.TLR16h <- res.TLR16h$log2FoldChange
names(gene_list.TLR16h) <- res.TLR16h$ENSG
  
set.seed(4735)
gseGO.TLR16h <- gseGO(gene_list.TLR16h,
                     OrgDb = org.Hs.eg.db,
                     keyType = "ENSEMBL",
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = TRUE)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseGO.TLR16h <- setReadable(gseGO.TLR16h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLR16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      BP 
## #...@keytype      ENSEMBL 
## #...@geneList     Named num [1:32287] 1.7 1.63 1.55 1.54 1.46 ...
##  - attr(*, "names")= chr [1:32287] "ENSG00000140297" "ENSG00000108448" "ENSG00000234719" "ENSG00000087842" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...401 enriched terms found
## 'data.frame':    401 obs. of  11 variables:
##  $ ID             : chr  "GO:0007059" "GO:0098813" "GO:0000280" "GO:0048285" ...
##  $ Description    : chr  "chromosome segregation" "nuclear chromosome segregation" "nuclear division" "organelle fission" ...
##  $ setSize        : int  404 297 413 459 417 220 435 477 181 322 ...
##  $ enrichmentScore: num  -0.761 -0.78 -0.733 -0.712 -0.714 ...
##  $ NES            : num  -2.32 -2.33 -2.25 -2.19 -2.19 ...
##  $ pvalue         : num  7.68e-35 2.55e-29 7.75e-29 1.62e-28 2.86e-26 ...
##  $ p.adjust       : num  3.70e-31 6.13e-26 1.24e-25 1.95e-25 2.76e-23 ...
##  $ qvalue         : num  3.15e-31 5.22e-26 1.06e-25 1.66e-25 2.35e-23 ...
##  $ rank           : num  2472 2572 2779 3019 3039 ...
##  $ leading_edge   : chr  "tags=47%, list=8%, signal=44%" "tags=47%, list=8%, signal=44%" "tags=41%, list=9%, signal=38%" "tags=40%, list=9%, signal=37%" ...
##  $ core_enrichment: chr  "KIF2A/KMT5A/GOLGA8A/MND1/CEP63/ACTL6A/ARID1A/SMARCC2/MSH5/CHMP2B/DDB1/SLC25A5/HAUS2/DIS3L2/PUM2/NUP62/NR3C1/RB1"| __truncated__ "ACTR3/BAG6/KIF2A/KMT5A/MND1/ACTL6A/ARID1A/SMARCC2/MSH5/CHMP2B/DDB1/DIS3L2/NUP62/RB1/SMARCE1/CDK5RAP2/SEH1L/HNRN"| __truncated__ "MLH1/MAPRE1/RAD1/CDC42/EPS8/CHEK1/ACTR3/BAG6/KIF2A/KMT5A/MND1/MSH5/CHMP2B/DDB1/RAD51D/DIS3L2/NUP62/RB1/CKS2/REE"| __truncated__ "TESMIN/RAB11A/REEP3/TENT4A/RALBP1/GDAP1/ARHGEF10/MLH1/MAPRE1/RAD1/CDC42/EPS8/CORO1C/CHEK1/ACTR3/BAG6/MYO19/KIF2"| __truncated__ ...
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseGO.TLR16h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)

TLR16h.go.plot <- nrow(gseGO.TLR16h) > 0
dotplot(gseGO.TLR16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseGO.TLR16h, showCategory = 5, color.params = list(foldChange = gene_list.TLR16h),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Warning: ggrepel: 3 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

gseGO.TLR16h.pwts <- pairwise_termsim(gseGO.TLR16h)

emapplot(gseGO.TLR16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

KEGG

ids <- bitr(names(gene_list.TLR16h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
## 'select()' returned 1:many mapping between keys and columns
ids <- ids[!duplicated(ids$ENSEMBL),]

gene_list.TLR16h.entrez <- gene_list.TLR16h[names(gene_list.TLR16h) %in% ids$ENSEMBL]
names(gene_list.TLR16h.entrez) <- ids$ENTREZID
gene_list.TLR16h.entrez <- gene_list.TLR16h.entrez[!duplicated(names(gene_list.TLR16h.entrez))]
set.seed(5923)
gseKEGG.TLR16h <- gseKEGG(gene_list.TLR16h.entrez,
                         organism = "hsa",
                         keyType = "ncbi-geneid",
                         eps = 0,
                         maxGSSize = 500,
                         minGSSize = 15,
                         seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseKEGG.TLR16h <- setReadable(gseKEGG.TLR16h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLR16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     hsa 
## #...@setType      KEGG 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 1.7 1.55 1.54 1.46 1.42 ...
##  - attr(*, "names")= chr [1:22414] "9245" "729978" "8544" "2069" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...42 enriched terms found
## 'data.frame':    42 obs. of  11 variables:
##  $ ID             : chr  "hsa04110" "hsa03010" "hsa04510" "hsa03030" ...
##  $ Description    : chr  "Cell cycle" "Ribosome" "Focal adhesion" "DNA replication" ...
##  $ setSize        : int  155 129 195 36 41 54 173 212 106 38 ...
##  $ enrichmentScore: num  -0.753 0.664 -0.621 -0.857 -0.843 ...
##  $ NES            : num  -2.27 2.12 -1.92 -2.1 -2.1 ...
##  $ pvalue         : num  4.04e-16 1.18e-09 6.94e-09 6.23e-08 6.54e-08 ...
##  $ p.adjust       : num  1.35e-13 1.98e-07 7.72e-07 3.12e-06 3.12e-06 ...
##  $ qvalue         : num  1.06e-13 1.56e-07 6.08e-07 2.46e-06 2.46e-06 ...
##  $ rank           : num  1963 6043 3758 2048 2101 ...
##  $ leading_edge   : chr  "tags=47%, list=9%, signal=43%" "tags=71%, list=27%, signal=52%" "tags=42%, list=17%, signal=35%" "tags=78%, list=9%, signal=71%" ...
##  $ core_enrichment: chr  "RB1/CCND1/TGFB2/CDC25B/EP300/ANAPC15/PDS5B/STAG1/ANAPC5/CCNH/PTTG1/NIPBL/ATM/CDT1/CDKN2D/BUB3/CREBBP/WEE1/ATR/C"| __truncated__ "MRPL23/RPS27/MRPL34/MRPL2/MRPL10/RPS15A/RPS19/RSL24D1/MRPS17/RPL28/RPL39/MRPL9/RPS28/MRPL16/RPL13A/MRPL36/MRPL1"| __truncated__ "IGF1R/PDGFB/COL4A1/ITGB5/PAK2/ITGB3/LAMA5/COL6A6/EMP2/LAMA4/LAMB1/ITGA5/CAV2/GSK3B/PIP5K1A/ITGA10/MYL12A/TNXB/M"| __truncated__ "RPA3/RNASEH1/POLE3/RFC4/RNASEH2A/POLD3/POLD1/RFC1/POLA2/RFC5/PCNA/PRIM1/RPA1/PRIM2/POLE2/MCM6/RFC2/FEN1/MCM4/MC"| __truncated__ ...
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseKEGG.TLR16h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)

TLR16h.kegg.plot <- nrow(gseKEGG.TLR16h) > 0
dotplot(gseKEGG.TLR16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseKEGG.TLR16h, showCategory = 5, color.params = list(foldChange = gene_list.TLR16h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseKEGG.TLR16h.pwts <- pairwise_termsim(gseKEGG.TLR16h)

emapplot(gseKEGG.TLR16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

pathview(gene.data = gene_list.TLR16h.entrez, pathway.id = gseKEGG.TLR16h@result$ID[1], species = "hsa")
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq
## Info: Writing image file hsa04110.pathview.png
filename <- paste0(gseKEGG.TLR16h@result$ID[1], ".pathview.png")

file.rename(filename, paste0(out.dir, filename))
## [1] TRUE
knitr::include_graphics(paste0(out.dir, filename))

Disease ontology

set.seed(6246)
gseDO.TLR16h <- gseDO(gene_list.TLR16h.entrez,
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseDO.TLR16h <- setReadable(gseDO.TLR16h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLR16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      DO 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 1.7 1.55 1.54 1.46 1.42 ...
##  - attr(*, "names")= chr [1:22414] "9245" "729978" "8544" "2069" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...96 enriched terms found
## 'data.frame':    96 obs. of  11 variables:
##  $ ID             : chr  "DOID:0080015" "DOID:5683" "DOID:2490" "DOID:1793" ...
##  $ Description    : chr  "physical disorder" "hereditary breast ovarian cancer syndrome" "congenital nervous system abnormality" "pancreatic cancer" ...
##  $ setSize        : int  302 213 135 453 460 264 255 144 144 29 ...
##  $ enrichmentScore: num  -0.581 -0.617 -0.665 -0.493 -0.489 ...
##  $ NES            : num  -1.87 -1.91 -1.97 -1.63 -1.62 ...
##  $ pvalue         : num  7.25e-10 5.94e-09 1.94e-08 3.13e-07 2.54e-07 ...
##  $ p.adjust       : num  6.96e-07 2.85e-06 6.19e-06 6.00e-05 6.00e-05 ...
##  $ qvalue         : num  5.29e-07 2.17e-06 4.71e-06 4.57e-05 4.57e-05 ...
##  $ rank           : num  3211 2757 3211 2639 2222 ...
##  $ leading_edge   : chr  "tags=33%, list=14%, signal=29%" "tags=40%, list=12%, signal=36%" "tags=40%, list=14%, signal=34%" "tags=26%, list=12%, signal=24%" ...
##  $ core_enrichment: chr  "NEDD4L/NOG/LAMB1/HMGA2/ARFGEF2/GOLGB1/STAG2/FGFR1/VANGL1/NUAK2/MAD1L1/MEST/CEP290/NSD1/DAW1/MECP2/TWSG1/RPGRIP1"| __truncated__ "CCNT1/NOLC1/EMSY/NPM1/MLH1/CHD8/CHEK1/GTF2I/DHX9/SP1/FASN/PALB2/NCOA2/RAD51D/HNRNPA2B1/COL1A1/RB1/CCND1/ERCC5/R"| __truncated__ "NEDD4L/LAMB1/ARFGEF2/STAG2/NUAK2/MAD1L1/NSD1/MECP2/TWSG1/RPGRIP1L/PNKP/CHEK1/DAG1/MTRR/MAP1B/CDK5RAP2/AUTS2/REL"| __truncated__ "RALBP1/PRKD1/NRP2/FANCC/MLH1/IKBKB/TSC2/LIG4/CDC42/BRAF/CHEK1/MTRR/APAF1/ILK/MAPK9/DNMT3B/ANXA2/CYCS/CTNNA1/GDI"| __truncated__ ...
## #...Citation
##   Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
##   R/Bioconductor package for Disease Ontology Semantic and Enrichment
##   analysis. Bioinformatics 2015, 31(4):608-609
write.csv(data.frame(gseDO.TLR16h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)

TLR16h.DO.plot <- nrow(gseDO.TLR16h) > 0
dotplot(gseDO.TLR16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseDO.TLR16h, showCategory = 5, color.params = list(foldChange = gene_list.TLR16h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Warning: ggrepel: 13 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

gseDO.TLR16h.pwts <- pairwise_termsim(gseDO.TLR16h)

emapplot(gseDO.TLR16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

TLR10 light vs. wildtype light 3h

Gene ontology

out.dir <- "./results/FA/TLR10_light_3h_vs_wildtype_light_3h/"

dir.create(out.dir)

res.TLRlightWT3h <- read.csv("results/DESeq/TLR10_light_3h_vs_wildtype_light_3h_unfiltered.csv") %>%
  as_tibble() %>%
  arrange(desc(log2FoldChange))

gene_list.TLRlightWT3h <- res.TLRlightWT3h$log2FoldChange
names(gene_list.TLRlightWT3h) <- res.TLRlightWT3h$ENSG
  
set.seed(4735)
gseGO.TLRlightWT3h <- gseGO(gene_list.TLRlightWT3h,
                     OrgDb = org.Hs.eg.db,
                     keyType = "ENSEMBL",
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = TRUE)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseGO.TLRlightWT3h <- setReadable(gseGO.TLRlightWT3h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLRlightWT3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      BP 
## #...@keytype      ENSEMBL 
## #...@geneList     Named num [1:32287] 9.67 3.47 2.11 1.99 1.96 ...
##  - attr(*, "names")= chr [1:32287] "ENSG00000174123" "ENSG00000164185" "ENSG00000123453" "ENSG00000225756" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...12 enriched terms found
## 'data.frame':    12 obs. of  11 variables:
##  $ ID             : chr  "GO:0030865" "GO:0002224" "GO:0030866" "GO:0032528" ...
##  $ Description    : chr  "cortical cytoskeleton organization" "toll-like receptor signaling pathway" "cortical actin cytoskeleton organization" "microvillus organization" ...
##  $ setSize        : int  50 63 42 24 58 213 83 65 56 68 ...
##  $ enrichmentScore: num  -0.997 0.998 -0.998 -0.999 -0.996 ...
##  $ NES            : num  -2.27 2.16 -2.22 -2.08 -2.29 ...
##  $ pvalue         : num  1.15e-06 1.63e-06 4.89e-06 5.39e-06 1.88e-05 ...
##  $ p.adjust       : num  0.00393 0.00393 0.00649 0.00649 0.01634 ...
##  $ qvalue         : num  0.0037 0.0037 0.00612 0.00612 0.0154 ...
##  $ rank           : num  2 7 2 10 2 2 2 6 6 6 ...
##  $ leading_edge   : chr  "tags=4%, list=0%, signal=4%" "tags=3%, list=0%, signal=3%" "tags=5%, list=0%, signal=5%" "tags=12%, list=0%, signal=13%" ...
##  $ core_enrichment: chr  "PLEK/LCP1" "TLR10/CD40" "PLEK/LCP1" "GLDN/TNIK/STK26" ...
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseGO.TLRlightWT3h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)

TLRlightWT3h.go.plot <- nrow(gseGO.TLRlightWT3h) > 0
dotplot(gseGO.TLRlightWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseGO.TLRlightWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT3h),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseGO.TLRlightWT3h.pwts <- pairwise_termsim(gseGO.TLRlightWT3h)

emapplot(gseGO.TLRlightWT3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

KEGG

ids <- bitr(names(gene_list.TLRlightWT3h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
## 'select()' returned 1:many mapping between keys and columns
ids <- ids[!duplicated(ids$ENSEMBL),]

gene_list.TLRlightWT3h.entrez <- gene_list.TLRlightWT3h[names(gene_list.TLRlightWT3h) %in% ids$ENSEMBL]
names(gene_list.TLRlightWT3h.entrez) <- ids$ENTREZID
gene_list.TLRlightWT3h.entrez <- gene_list.TLRlightWT3h.entrez[!duplicated(names(gene_list.TLRlightWT3h.entrez))]
set.seed(5923)
gseKEGG.TLRlightWT3h <- gseKEGG(gene_list.TLRlightWT3h.entrez,
                         organism = "hsa",
                         keyType = "ncbi-geneid",
                         eps = 0,
                         maxGSSize = 500,
                         minGSSize = 15,
                         seed = T)
## preparing geneSet collections...
## GSEA analysis...
## no term enriched under specific pvalueCutoff...
gseKEGG.TLRlightWT3h <- setReadable(gseKEGG.TLRlightWT3h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLRlightWT3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     hsa 
## #...@setType      KEGG 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 9.67 3.47 2.11 1.99 1.96 ...
##  - attr(*, "names")= chr [1:22414] "81793" "133923" "1757" "138948" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...0 enriched terms found
## 'data.frame':    0 obs. of  8 variables:
##  $ ID             : chr 
##  $ Description    : chr 
##  $ setSize        : int 
##  $ enrichmentScore: num 
##  $ NES            : num 
##  $ pvalue         : num 
##  $ p.adjust       : num 
##  $ qvalue         : num 
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseKEGG.TLRlightWT3h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)

TLRlightWT3h.kegg.plot <- nrow(gseKEGG.TLRlightWT3h) > 0
dotplot(gseKEGG.TLRlightWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
cnetplot(gseKEGG.TLRlightWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT3h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
gseKEGG.TLRlightWT3h.pwts <- pairwise_termsim(gseKEGG.TLRlightWT3h)

emapplot(gseKEGG.TLRlightWT3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
pathview(gene.data = gene_list.TLRlightWT3h.entrez, pathway.id = gseKEGG.TLRlightWT3h@result$ID[1], species = "hsa")

filename <- paste0(gseKEGG.TLRlightWT3h@result$ID[1], ".pathview.png")

file.rename(filename, paste0(out.dir, filename))

knitr::include_graphics(paste0(out.dir, filename))

Disease ontology

set.seed(6246)
gseDO.TLRlightWT3h <- gseDO(gene_list.TLRlightWT3h.entrez,
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = T)
## preparing geneSet collections...
## GSEA analysis...
## no term enriched under specific pvalueCutoff...
gseDO.TLRlightWT3h <- setReadable(gseDO.TLRlightWT3h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLRlightWT3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      DO 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 9.67 3.47 2.11 1.99 1.96 ...
##  - attr(*, "names")= chr [1:22414] "81793" "133923" "1757" "138948" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...0 enriched terms found
## 'data.frame':    0 obs. of  8 variables:
##  $ ID             : chr 
##  $ Description    : chr 
##  $ setSize        : int 
##  $ enrichmentScore: num 
##  $ NES            : num 
##  $ pvalue         : num 
##  $ p.adjust       : num 
##  $ qvalue         : num 
## #...Citation
##   Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
##   R/Bioconductor package for Disease Ontology Semantic and Enrichment
##   analysis. Bioinformatics 2015, 31(4):608-609
write.csv(data.frame(gseDO.TLRlightWT3h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)

TLRlightWT3h.DO.plot <- nrow(gseDO.TLRlightWT3h) > 0
dotplot(gseDO.TLRlightWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
cnetplot(gseDO.TLRlightWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT3h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
gseDO.TLRlightWT3h.pwts <- pairwise_termsim(gseDO.TLRlightWT3h)

emapplot(gseDO.TLRlightWT3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")

TLR10 light vs. wildtype light 16h

Gene ontology

out.dir <- "./results/FA/TLR10_light_16h_vs_wildtype_light_16h/"

dir.create(out.dir)

res.TLRlightWT16h <- read.csv("results/DESeq/TLR10_light_16h_vs_wildtype_light_16h_unfiltered.csv") %>%
  as_tibble() %>%
  arrange(desc(log2FoldChange))

gene_list.TLRlightWT16h <- res.TLRlightWT16h$log2FoldChange
names(gene_list.TLRlightWT16h) <- res.TLRlightWT16h$ENSG
  
set.seed(6548)
gseGO.TLRlightWT16h <- gseGO(gene_list.TLRlightWT16h,
                     OrgDb = org.Hs.eg.db,
                     keyType = "ENSEMBL",
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = TRUE)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseGO.TLRlightWT16h <- setReadable(gseGO.TLRlightWT16h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLRlightWT16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      BP 
## #...@keytype      ENSEMBL 
## #...@geneList     Named num [1:32287] 6.72 2.57 2.19 1.92 1.57 ...
##  - attr(*, "names")= chr [1:32287] "ENSG00000174123" "ENSG00000225756" "ENSG00000139292" "ENSG00000123453" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...3 enriched terms found
## 'data.frame':    3 obs. of  11 variables:
##  $ ID             : chr  "GO:2000052" "GO:2000095" "GO:0072111"
##  $ Description    : chr  "positive regulation of non-canonical Wnt signaling pathway" "regulation of Wnt signaling pathway, planar cell polarity pathway" "cell proliferation involved in kidney development"
##  $ setSize        : int  15 16 21
##  $ enrichmentScore: num  -0.986 -0.986 -0.987
##  $ NES            : num  -1.87 -1.91 -2.01
##  $ pvalue         : num  1.66e-06 4.25e-06 1.12e-05
##  $ p.adjust       : num  0.00798 0.01024 0.01792
##  $ qvalue         : num  0.00678 0.0087 0.01523
##  $ rank           : num  2 148 2
##  $ leading_edge   : chr  "tags=13%, list=0%, signal=13%" "tags=12%, list=0%, signal=12%" "tags=10%, list=0%, signal=10%"
##  $ core_enrichment: chr  "WNT5A/GPC3" "SFRP2/GPC3" "ITGB3/GPC3"
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseGO.TLRlightWT16h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)

TLRlightWT16h.go.plot <- nrow(gseGO.TLRlightWT16h) > 0
dotplot(gseGO.TLRlightWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseGO.TLRlightWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT16h),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseGO.TLRlightWT16h.pwts <- pairwise_termsim(gseGO.TLRlightWT16h)

emapplot(gseGO.TLRlightWT16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

KEGG

ids <- bitr(names(gene_list.TLRlightWT16h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
## 'select()' returned 1:many mapping between keys and columns
ids <- ids[!duplicated(ids$ENSEMBL),]

gene_list.TLRlightWT16h.entrez <- gene_list.TLRlightWT16h[names(gene_list.TLRlightWT16h) %in% ids$ENSEMBL]
names(gene_list.TLRlightWT16h.entrez) <- ids$ENTREZID
gene_list.TLRlightWT16h.entrez <- gene_list.TLRlightWT16h.entrez[!duplicated(names(gene_list.TLRlightWT16h.entrez))]
set.seed(5923)
gseKEGG.TLRlightWT16h <- gseKEGG(gene_list.TLRlightWT16h.entrez,
                         organism = "hsa",
                         keyType = "ncbi-geneid",
                         eps = 0,
                         maxGSSize = 500,
                         minGSSize = 15,
                         seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseKEGG.TLRlightWT16h <- setReadable(gseKEGG.TLRlightWT16h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLRlightWT16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     hsa 
## #...@setType      KEGG 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 6.72 2.57 2.19 1.92 1.57 ...
##  - attr(*, "names")= chr [1:22414] "81793" "138948" "8549" "1757" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...1 enriched terms found
## 'data.frame':    1 obs. of  11 variables:
##  $ ID             : chr "hsa03010"
##  $ Description    : chr "Ribosome"
##  $ setSize        : int 129
##  $ enrichmentScore: num 0.746
##  $ NES            : num 2.25
##  $ pvalue         : num 9.41e-06
##  $ p.adjust       : num 0.00314
##  $ qvalue         : num 0.00276
##  $ rank           : int 3268
##  $ leading_edge   : chr "tags=69%, list=15%, signal=59%"
##  $ core_enrichment: chr "MRPS17/RPS19/MRPL34/RPS15A/RPS29/RPLP1/RPL28/UBA52/RPL13A/RPS27L/RPL13/MRPS21/RPL37A/MRPL36/RPL27A/RPLP2/MRPL4/"| __truncated__
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseKEGG.TLRlightWT16h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)

TLRlightWT16h.kegg.plot <- nrow(gseKEGG.TLRlightWT16h) > 0
dotplot(gseKEGG.TLRlightWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseKEGG.TLRlightWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT16h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseKEGG.TLRlightWT16h.pwts <- pairwise_termsim(gseKEGG.TLRlightWT16h)

emapplot(gseKEGG.TLRlightWT16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")

pathview(gene.data = gene_list.TLRlightWT16h.entrez, pathway.id = gseKEGG.TLRlightWT16h@result$ID[1], species = "hsa")
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq
## Info: Writing image file hsa03010.pathview.png
filename <- paste0(gseKEGG.TLRlightWT16h@result$ID[1], ".pathview.png")

file.rename(filename, paste0(out.dir, filename))
## [1] TRUE
knitr::include_graphics(paste0(out.dir, filename))

Disease ontology

set.seed(6246)
gseDO.TLRlightWT16h <- gseDO(gene_list.TLRlightWT16h.entrez,
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseDO.TLRlightWT16h <- setReadable(gseDO.TLRlightWT16h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLRlightWT16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      DO 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 6.72 2.57 2.19 1.92 1.57 ...
##  - attr(*, "names")= chr [1:22414] "81793" "138948" "8549" "1757" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...33 enriched terms found
## 'data.frame':    33 obs. of  11 variables:
##  $ ID             : chr  "DOID:3113" "DOID:4481" "DOID:2021" "DOID:3594" ...
##  $ Description    : chr  "papillary carcinoma" "allergic rhinitis" "placenta cancer" "choriocarcinoma" ...
##  $ setSize        : int  44 87 33 33 154 118 51 63 33 128 ...
##  $ enrichmentScore: num  -0.962 0.813 -0.973 -0.973 0.706 ...
##  $ NES            : num  -2.41 2.36 -2.32 -2.32 2.2 ...
##  $ pvalue         : num  1.38e-05 1.14e-05 7.06e-06 7.06e-06 1.07e-05 ...
##  $ p.adjust       : num  0.00264 0.00264 0.00264 0.00264 0.00264 ...
##  $ qvalue         : num  0.00237 0.00237 0.00237 0.00237 0.00237 ...
##  $ rank           : num  2 5 2 2 5 5 333 2 2 5 ...
##  $ leading_edge   : chr  "tags=5%, list=0%, signal=5%" "tags=2%, list=0%, signal=2%" "tags=6%, list=0%, signal=6%" "tags=6%, list=0%, signal=6%" ...
##  $ core_enrichment: chr  "MET/GPC3" "TLR10/CD40" "MTOR/GPC3" "MTOR/GPC3" ...
## #...Citation
##   Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
##   R/Bioconductor package for Disease Ontology Semantic and Enrichment
##   analysis. Bioinformatics 2015, 31(4):608-609
write.csv(data.frame(gseDO.TLRlightWT16h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)

TLRlightWT16h.DO.plot <- nrow(gseDO.TLRlightWT16h) > 0
dotplot(gseDO.TLRlightWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseDO.TLRlightWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT16h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseDO.TLRlightWT16h.pwts <- pairwise_termsim(gseDO.TLRlightWT16h)

emapplot(gseDO.TLRlightWT16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

TLR10 dark vs. wildtype light 3h

Gene ontology

out.dir <- "./results/FA/TLR10_dark_3h_vs_wildtype_light_3h/"

dir.create(out.dir)

res.TLRdarkWT3h <- read.csv("results/DESeq/TLR10_dark_3h_vs_wildtype_light_3h_unfiltered.csv") %>%
  as_tibble() %>%
  arrange(desc(log2FoldChange))

gene_list.TLRdarkWT3h <- res.TLRdarkWT3h$log2FoldChange
names(gene_list.TLRdarkWT3h) <- res.TLRdarkWT3h$ENSG
  
set.seed(4735)
gseGO.TLRdarkWT3h <- gseGO(gene_list.TLRdarkWT3h,
                     OrgDb = org.Hs.eg.db,
                     keyType = "ENSEMBL",
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = TRUE)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseGO.TLRdarkWT3h <- setReadable(gseGO.TLRdarkWT3h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLRdarkWT3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      BP 
## #...@keytype      ENSEMBL 
## #...@geneList     Named num [1:32287] 28.9 27.5 19.4 11.2 10.1 ...
##  - attr(*, "names")= chr [1:32287] "ENSG00000157851" "ENSG00000132182" "ENSG00000165588" "ENSG00000130294" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...17 enriched terms found
## 'data.frame':    17 obs. of  11 variables:
##  $ ID             : chr  "GO:0006986" "GO:0031667" "GO:0009991" "GO:0070059" ...
##  $ Description    : chr  "response to unfolded protein" "response to nutrient levels" "response to extracellular stimulus" "intrinsic apoptotic signaling pathway in response to endoplasmic reticulum stress" ...
##  $ setSize        : int  134 437 464 61 200 202 365 495 296 170 ...
##  $ enrichmentScore: num  -0.764 -0.604 -0.599 -0.852 -0.702 ...
##  $ NES            : num  -2 -1.78 -1.77 -2.01 -1.92 ...
##  $ pvalue         : num  3.37e-06 1.25e-06 2.42e-06 1.18e-05 1.05e-05 ...
##  $ p.adjust       : num  0.00542 0.00542 0.00542 0.01139 0.01139 ...
##  $ qvalue         : num  0.00526 0.00526 0.00526 0.01106 0.01106 ...
##  $ rank           : num  1924 3563 3563 767 3196 ...
##  $ leading_edge   : chr  "tags=24%, list=6%, signal=23%" "tags=26%, list=11%, signal=24%" "tags=26%, list=11%, signal=24%" "tags=25%, list=2%, signal=24%" ...
##  $ core_enrichment: chr  "RHBDD1/DDRGK1/CREB3L4/ERN1/CDK5RAP3/HSPA6/PARP8/HSPA8/BFAR/OS9/CTH/EDEM2/ATF3/DNAJB2/CREBRF/RACK1/MANF/SERPINH1"| __truncated__ "SORL1/SPX/COMT/PRKAB2/TSSK6/MAP1LC3A/GABARAPL1/RXRA/LRRK2/VPS41/ULK1/MLST8/FES/APOE/PARP2/PPARA/PLD1/TSPO/ACE/L"| __truncated__ "SORL1/SPX/COMT/PRKAB2/TSSK6/MAP1LC3A/GABARAPL1/RXRA/NOCT/LRRK2/VPS41/ULK1/MLST8/FES/APOE/PARP2/PPARA/PLD1/TSPO/"| __truncated__ "HYOU1/CASP4/SYVN1/BCL2/BBC3/ERP29/CEBPB/TRIB3/PMAIP1/XBP1/ATF4/HSPA1A/HERPUD1/CHAC1/DDIT3" ...
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseGO.TLRdarkWT3h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)

TLRdarkWT3h.go.plot <- nrow(gseGO.TLRdarkWT3h) > 0
dotplot(gseGO.TLRdarkWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseGO.TLRdarkWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT3h),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseGO.TLRdarkWT3h.pwts <- pairwise_termsim(gseGO.TLRdarkWT3h)

emapplot(gseGO.TLRdarkWT3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

KEGG

ids <- bitr(names(gene_list.TLRdarkWT3h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
## 'select()' returned 1:many mapping between keys and columns
ids <- ids[!duplicated(ids$ENSEMBL),]

gene_list.TLRdarkWT3h.entrez <- gene_list.TLRdarkWT3h[names(gene_list.TLRdarkWT3h) %in% ids$ENSEMBL]
names(gene_list.TLRdarkWT3h.entrez) <- ids$ENTREZID
gene_list.TLRdarkWT3h.entrez <- gene_list.TLRdarkWT3h.entrez[!duplicated(names(gene_list.TLRdarkWT3h.entrez))]
set.seed(5923)
gseKEGG.TLRdarkWT3h <- gseKEGG(gene_list.TLRdarkWT3h.entrez,
                         organism = "hsa",
                         keyType = "ncbi-geneid",
                         eps = 0,
                         maxGSSize = 500,
                         minGSSize = 15,
                         seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseKEGG.TLRdarkWT3h <- setReadable(gseKEGG.TLRdarkWT3h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLRdarkWT3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     hsa 
## #...@setType      KEGG 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 28.9 27.5 19.4 11.2 10.1 ...
##  - attr(*, "names")= chr [1:22414] "56896" "23225" "5015" "547" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...3 enriched terms found
## 'data.frame':    3 obs. of  11 variables:
##  $ ID             : chr  "hsa03010" "hsa04141" "hsa04612"
##  $ Description    : chr  "Ribosome" "Protein processing in endoplasmic reticulum" "Antigen processing and presentation"
##  $ setSize        : int  129 164 57
##  $ enrichmentScore: num  -0.735 -0.696 -0.812
##  $ NES            : num  -2.05 -2.01 -2.02
##  $ pvalue         : num  6.10e-06 3.56e-06 1.33e-04
##  $ p.adjust       : num  0.00102 0.00102 0.01475
##  $ qvalue         : num  0.000979 0.000979 0.014182
##  $ rank           : num  2463 2245 1726
##  $ leading_edge   : chr  "tags=59%, list=11%, signal=53%" "tags=29%, list=10%, signal=27%" "tags=32%, list=8%, signal=29%"
##  $ core_enrichment: chr  "RPS15A/RPL36AL/RPL41/MRPS15/RPL24/RPL35/RPL39/RPL23A/MRPL10/RPS5/RPL10/RPL6/RPLP0/RPS27A/RPL15/RPS7/RPL22/RPL18"| __truncated__ "SSR2/DNAJC3/DNAJC1/TRAF2/UBXN6/TUSC3/LMAN1/ERO1B/MAN1B1/TXNDC5/PDIA4/BAG2/SIL1/ERN1/HSPA6/PRKCSH/CRYAB/DNAJB11/"| __truncated__ "NFYA/HLA-E/HLA-F/HSPA6/HSPA8/HLA-B/NFYB/B2M/CTSL/HLA-C/PSME1/HLA-A/CREB1/TAP1/HSPA5/HSPA1L/HSPA1A/HSPA1B"
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseKEGG.TLRdarkWT3h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)

TLRdarkWT3h.kegg.plot <- nrow(gseKEGG.TLRdarkWT3h) > 0
dotplot(gseKEGG.TLRdarkWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseKEGG.TLRdarkWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT3h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseKEGG.TLRdarkWT3h.pwts <- pairwise_termsim(gseKEGG.TLRdarkWT3h)

emapplot(gseKEGG.TLRdarkWT3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

pathview(gene.data = gene_list.TLRdarkWT3h.entrez, pathway.id = gseKEGG.TLRdarkWT3h@result$ID[1], species = "hsa")
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq
## Info: Writing image file hsa03010.pathview.png
filename <- paste0(gseKEGG.TLRdarkWT3h@result$ID[1], ".pathview.png")

file.rename(filename, paste0(out.dir, filename))
## [1] TRUE
knitr::include_graphics(paste0(out.dir, filename))

Disease ontology

set.seed(6246)
gseDO.TLRdarkWT3h <- gseDO(gene_list.TLRdarkWT3h.entrez,
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseDO.TLRdarkWT3h <- setReadable(gseDO.TLRdarkWT3h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLRdarkWT3h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      DO 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 28.9 27.5 19.4 11.2 10.1 ...
##  - attr(*, "names")= chr [1:22414] "56896" "23225" "5015" "547" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...8 enriched terms found
## 'data.frame':    8 obs. of  11 variables:
##  $ ID             : chr  "DOID:11949" "DOID:12306" "DOID:3021" "DOID:10825" ...
##  $ Description    : chr  "Creutzfeldt-Jakob disease" "vitiligo" "acute kidney failure" "essential hypertension" ...
##  $ setSize        : int  23 70 99 111 289 459 79 470
##  $ enrichmentScore: num  -0.915 -0.783 -0.725 -0.711 -0.554 ...
##  $ NES            : num  -2.01 -2 -1.93 -1.92 -1.68 ...
##  $ pvalue         : num  1.50e-04 1.93e-04 8.59e-05 1.31e-04 2.01e-04 ...
##  $ p.adjust       : num  0.0321 0.0321 0.0321 0.0321 0.0321 ...
##  $ qvalue         : num  0.0285 0.0285 0.0285 0.0285 0.0285 ...
##  $ rank           : num  243 2235 1728 2636 1861 ...
##  $ leading_edge   : chr  "tags=13%, list=1%, signal=13%" "tags=24%, list=10%, signal=22%" "tags=22%, list=8%, signal=21%" "tags=23%, list=12%, signal=20%" ...
##  $ core_enrichment: chr  "ABCB1/SNCA/PTGS2" "HGF/IL6/FAS/TXNDC5/GSTP1/PSMB8/HLA-B/PSMB9/HLA-C/CAT/HLA-A/TAP1/TSLP/NFE2L2/XBP1/PTGS2/HSPA1A" "PPARG/GFER/LGALS3/FIS1/HMOX1/HSPA8/B2M/HYOU1/CLU/CASP1/SLC22A1/PTAFR/TNFRSF1B/CDKN1B/NQO1/BCL2/IL6R/NFE2L2/ARG2"| __truncated__ "FMO3/ADA/MMP2/GSTM3/ERAP1/TGFB3/HGF/IL6/ATP5PF/CST3/PPARG/HMOX1/MTHFR/KYNU/CTH/PSMB9/GDF15/ALDH2/CAT/CLU/RGS5/M"| __truncated__ ...
## #...Citation
##   Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
##   R/Bioconductor package for Disease Ontology Semantic and Enrichment
##   analysis. Bioinformatics 2015, 31(4):608-609
write.csv(data.frame(gseDO.TLRdarkWT3h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)

TLRdarkWT3h.DO.plot <- nrow(gseDO.TLRdarkWT3h) > 0
dotplot(gseDO.TLRdarkWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseDO.TLRdarkWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT3h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseDO.TLRdarkWT3h.pwts <- pairwise_termsim(gseDO.TLRdarkWT3h)

emapplot(gseDO.TLRdarkWT3h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

TLR10 dark vs. wildtype light 16h

Gene ontology

out.dir <- "./results/FA/TLR10_dark_16h_vs_wildtype_light_16h/"

dir.create(out.dir)

res.TLRdarkWT16h <- read.csv("results/DESeq/TLR10_dark_16h_vs_wildtype_light_16h_unfiltered.csv") %>%
  as_tibble() %>%
  arrange(desc(log2FoldChange))

gene_list.TLRdarkWT16h <- res.TLRdarkWT16h$log2FoldChange
names(gene_list.TLRdarkWT16h) <- res.TLRdarkWT16h$ENSG
  
set.seed(4735)
gseGO.TLRdarkWT16h <- gseGO(gene_list.TLRdarkWT16h,
                     OrgDb = org.Hs.eg.db,
                     keyType = "ENSEMBL",
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = TRUE)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseGO.TLRdarkWT16h <- setReadable(gseGO.TLRdarkWT16h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLRdarkWT16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      BP 
## #...@keytype      ENSEMBL 
## #...@geneList     Named num [1:32287] 7.73 3.18 2.8 2.75 2.74 ...
##  - attr(*, "names")= chr [1:32287] "ENSG00000174123" "ENSG00000162496" "ENSG00000139292" "ENSG00000225756" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...32 enriched terms found
## 'data.frame':    32 obs. of  11 variables:
##  $ ID             : chr  "GO:0007059" "GO:0098813" "GO:0051983" "GO:0000819" ...
##  $ Description    : chr  "chromosome segregation" "nuclear chromosome segregation" "regulation of chromosome segregation" "sister chromatid segregation" ...
##  $ setSize        : int  404 297 125 220 53 72 61 58 47 47 ...
##  $ enrichmentScore: num  0.629 0.656 0.757 0.679 0.837 ...
##  $ NES            : num  1.77 1.81 1.89 1.8 1.9 ...
##  $ pvalue         : num  4.00e-06 8.32e-06 1.29e-05 1.86e-05 4.33e-05 ...
##  $ p.adjust       : num  0.0193 0.02 0.0208 0.0224 0.0232 ...
##  $ qvalue         : num  0.0184 0.0191 0.0198 0.0214 0.0221 ...
##  $ rank           : num  3353 4105 3299 3995 3256 ...
##  $ leading_edge   : chr  "tags=44%, list=10%, signal=40%" "tags=49%, list=13%, signal=43%" "tags=53%, list=10%, signal=48%" "tags=51%, list=12%, signal=45%" ...
##  $ core_enrichment: chr  "TTN/FBXO5/CENPK/DPF3/AURKB/NEK7/BRIP1/RMDN1/GOLGA8A/PRICKLE1/CENPI/DSN1/CHMP3/CENPQ/KNTC1/SKA2/RMI2/CENPU/KIF23"| __truncated__ "TTN/FBXO5/CENPK/DPF3/AURKB/BRIP1/RMDN1/PRICKLE1/CENPI/DSN1/CHMP3/CENPQ/KNTC1/RMI2/KIF23/PSMC3IP/CDC6/PMF1/CDK1/"| __truncated__ "FBXO5/DPF3/AURKB/KNTC1/RMI2/CDC6/CDK1/RAD18/HASPIN/NCAPG/SPDL1/BUB1B/GEN1/SMC2/NCAPH/ZWINT/ANAPC5/TTK/NCAPG2/SP"| __truncated__ "TTN/FBXO5/CENPK/DPF3/AURKB/RMDN1/PRICKLE1/CENPI/DSN1/CHMP3/KNTC1/RMI2/KIF23/CDC6/CDK1/RACGAP1/HASPIN/NCAPG/SPDL"| __truncated__ ...
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseGO.TLRdarkWT16h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)

TLRdarkWT16h.go.plot <- nrow(gseGO.TLRdarkWT16h) > 0
dotplot(gseGO.TLRdarkWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseGO.TLRdarkWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT16h),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseGO.TLRdarkWT16h.pwts <- pairwise_termsim(gseGO.TLRdarkWT16h)

emapplot(gseGO.TLRdarkWT16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

KEGG

ids <- bitr(names(gene_list.TLRdarkWT16h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
## 'select()' returned 1:many mapping between keys and columns
ids <- ids[!duplicated(ids$ENSEMBL),]

gene_list.TLRdarkWT16h.entrez <- gene_list.TLRdarkWT16h[names(gene_list.TLRdarkWT16h) %in% ids$ENSEMBL]
names(gene_list.TLRdarkWT16h.entrez) <- ids$ENTREZID
gene_list.TLRdarkWT16h.entrez <- gene_list.TLRdarkWT16h.entrez[!duplicated(names(gene_list.TLRdarkWT16h.entrez))]
set.seed(5923)
gseKEGG.TLRdarkWT16h <- gseKEGG(gene_list.TLRdarkWT16h.entrez,
                         organism = "hsa",
                         keyType = "ncbi-geneid",
                         eps = 0,
                         maxGSSize = 500,
                         minGSSize = 15,
                         seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseKEGG.TLRdarkWT16h <- setReadable(gseKEGG.TLRdarkWT16h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLRdarkWT16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     hsa 
## #...@setType      KEGG 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 7.73 3.18 2.8 2.75 2.74 ...
##  - attr(*, "names")= chr [1:22414] "81793" "9249" "8549" "138948" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...2 enriched terms found
## 'data.frame':    2 obs. of  11 variables:
##  $ ID             : chr  "hsa05204" "hsa00190"
##  $ Description    : chr  "Chemical carcinogenesis - DNA adducts" "Oxidative phosphorylation"
##  $ setSize        : int  38 131
##  $ enrichmentScore: num  -0.789 0.641
##  $ NES            : num  -1.89 1.74
##  $ pvalue         : num  0.000103 0.000285
##  $ p.adjust       : num  0.0346 0.0476
##  $ qvalue         : num  0.0316 0.0435
##  $ rank           : num  1887 3982
##  $ leading_edge   : chr  "tags=32%, list=8%, signal=29%" "tags=56%, list=18%, signal=46%"
##  $ core_enrichment: chr  "GSTM4/CBR1/CYP1A1/GSTA4/GSTM3/GSTM2/MGST1/AKR1C2/EPHX1/CYP1B1/HSD11B1/PTGS2" "ATP6V0E2/ATP5F1D/COX3/COX2/ATP5F1C/NDUFB8/ATP5PO/NDUFA10/ATP5MG/NDUFB1/ATP5F1B/NDUFA11/ATP5F1A/COX6C/COX7B/NDUF"| __truncated__
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
write.csv(data.frame(gseKEGG.TLRdarkWT16h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)

TLRdarkWT16h.kegg.plot <- nrow(gseKEGG.TLRdarkWT16h) > 0
dotplot(gseKEGG.TLRdarkWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseKEGG.TLRdarkWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT16h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseKEGG.TLRdarkWT16h.pwts <- pairwise_termsim(gseKEGG.TLRdarkWT16h)

emapplot(gseKEGG.TLRdarkWT16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

pathview(gene.data = gene_list.TLRdarkWT16h.entrez, pathway.id = gseKEGG.TLRdarkWT16h@result$ID[1], species = "hsa")
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq
## Info: Writing image file hsa05204.pathview.png
filename <- paste0(gseKEGG.TLRdarkWT16h@result$ID[1], ".pathview.png")

file.rename(filename, paste0(out.dir, filename))
## [1] TRUE
knitr::include_graphics(paste0(out.dir, filename))

Disease ontology

set.seed(6246)
gseDO.TLRdarkWT16h <- gseDO(gene_list.TLRdarkWT16h.entrez,
                     eps = 0,
                     maxGSSize = 500,
                     minGSSize = 15,
                     seed = T)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
gseDO.TLRdarkWT16h <- setReadable(gseDO.TLRdarkWT16h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLRdarkWT16h
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     Homo sapiens 
## #...@setType      DO 
## #...@keytype      ENTREZID 
## #...@geneList     Named num [1:22414] 7.73 3.18 2.8 2.75 2.74 ...
##  - attr(*, "names")= chr [1:22414] "81793" "9249" "8549" "138948" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...3 enriched terms found
## 'data.frame':    3 obs. of  11 variables:
##  $ ID             : chr  "DOID:9974" "DOID:2559" "DOID:974"
##  $ Description    : chr  "drug dependence" "opiate dependence" "upper respiratory tract disease"
##  $ setSize        : int  18 17 154
##  $ enrichmentScore: num  0.941 0.941 0.651
##  $ NES            : num  1.89 1.86 1.8
##  $ pvalue         : num  5.02e-06 1.57e-05 1.46e-04
##  $ p.adjust       : num  0.00482 0.00754 0.04667
##  $ qvalue         : num  0.00445 0.00697 0.04313
##  $ rank           : int  21 21 663
##  $ leading_edge   : chr  "tags=11%, list=0%, signal=11%" "tags=12%, list=0%, signal=12%" "tags=6%, list=3%, signal=6%"
##  $ core_enrichment: chr  "KCNJ6/ABAT" "KCNJ6/ABAT" "TLR10/KCNJ6/BDNF/EDN1/CD40/MCAM/NLRP3/TLR4/TGFBR2/ADM"
## #...Citation
##   Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
##   R/Bioconductor package for Disease Ontology Semantic and Enrichment
##   analysis. Bioinformatics 2015, 31(4):608-609
write.csv(data.frame(gseDO.TLRdarkWT16h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)

TLRdarkWT16h.DO.plot <- nrow(gseDO.TLRdarkWT16h) > 0
dotplot(gseDO.TLRdarkWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) + 
  facet_grid(~.sign) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cnetplot(gseDO.TLRdarkWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT16h.entrez),
         cex.params = list(category_label = 0.6, gene_label = 0.4)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

gseDO.TLRdarkWT16h.pwts <- pairwise_termsim(gseDO.TLRdarkWT16h)

emapplot(gseDO.TLRdarkWT16h.pwts, force = 1.5,
         cex.params = list(category_label = 0.8)) +
  scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.